function [phi] = simiImplicitmove2(phi,mphi,nelx,nely,xkesi,ykesi,xyita,yyita,J,dkesi,dyita,dt,kappa)
%SIMIIMPLICITMOVE Summary of this function goes here
%   Detailed explanation goes here
  
    p = (xkesi.^2+ykesi.^2)./J;
    q = (xyita.^2+yyita.^2)./J;
    r = (xkesi.*xyita+ykesi.*yyita)./J;
    
    ne = nelx * nely;
    A = sparse(ne);
    for i = 2 : nelx-1
        for j = 2: nely - 1
            A(node(i,j),node(i,j)) = 1 - kappa*dt/J(i,j)*(-1/(2*dkesi^2)*(q(i+1,j)-q(i-1,j))+1/(2*dkesi*dyita)*(r(i+1,j)-r(i-1,j))...
                +1/(2*dkesi*dyita)*(r(i,j+1)-r(i,j-1))-1/(2*dyita^2)*(p(i,j+1)-p(i,j-1)));
            A(node(i,j),node(i+1,j)) = -kappa*dt/J(i,j)*(1/(2*dkesi^2)*(q(i+1,j)+q(i,j))-1/(2*dkesi*dyita)*(r(i,j+1)+r(i,j)));
            A(node(i,j),node(i,j+1)) = -kappa*dt/J(i,j)*(1/(2*dyita^2)*(p(i,j+1)+p(i,j))-1/(2*dkesi*dyita)*(r(i+1,j)+r(i,j)));
            A(node(i,j),node(i-1,j)) = -kappa*dt/J(i,j)*(1/(2*dkesi^2)*(q(i,j)+q(i-1,j))-1/(2*dkesi*dyita)*(r(i,j)+r(i,j-1)));
            A(node(i,j),node(i,j-1)) = -kappa*dt/J(i,j)*(1/(2*dyita^2)*(p(i,j)+p(i,j-1))-1/(2*dkesi*dyita)*(r(i,j)+r(i-1,j)));        
        end
        A(node(i,1),node(i,1)) = 1 - kappa*dt/J(i,1)*(-1/(2*dkesi^2)*(q(i+1,1)-q(i-1,1))+1/(2*dkesi*dyita)*(r(i+1,1)-r(i-1,1))...
                +1/(2*dkesi*dyita)*(r(i,2)-r(i,1))-1/(2*dyita^2)*(p(i,2)-p(i,1)));
        A(node(i,1),node(i,2)) = -2*kappa*dt/J(i,1)*(1/(2*dyita^2)*(p(i,2)+p(i,1))-1/(2*dkesi*dyita)*(r(i+1,1)+r(i,1)));
        A(node(i,1),node(i+1,1)) = -kappa*dt/J(i,1)*(1/(2*dkesi^2)*(q(i+1,1)+q(i,1))-1/(2*dkesi*dyita)*(r(i,2)+r(i,1)));
        A(node(i,1),node(i-1,1)) = -kappa*dt/J(i,1)*(1/(2*dkesi^2)*(q(i,1)+q(i-1,1))-1/(2*dkesi*dyita)*(r(i,1)+r(i,1)));
    
        A(node(i,nely),node(i,nely)) = 1 - kappa*dt/J(i,nely)*(-1/(2*dkesi^2)*(q(i+1,nely)-q(i-1,nely))+1/(2*dkesi*dyita)*(r(i+1,nely)-r(i-1,nely))...
                +1/(2*dkesi*dyita)*(r(i,nely)-r(i,nely-1))-1/(2*dyita^2)*(p(i,nely)-p(i,nely-1)));
        A(node(i,nely),node(i,nely-1)) = -2*kappa*dt/J(i,nely)*(1/(2*dyita^2)*(p(i,nely)+p(i,nely-1))-1/(2*dkesi*dyita)*(r(i,nely)+r(i-1,nely)));
        A(node(i,nely),node(i+1,nely)) = -kappa*dt/J(i,nely)*(1/(2*dkesi^2)*(q(i+1,nely)+q(i,nely))-1/(2*dkesi*dyita)*(r(i,nely)+r(i,nely)));
        A(node(i,nely),node(i-1,nely)) = -kappa*dt/J(i,nely)*(1/(2*dkesi^2)*(q(i,nely)+q(i-1,nely))-1/(2*dkesi*dyita)*(r(i,nely)+r(i,nely-1)));
    end
    
    for j = 2 : nely - 1
    A(node(1,j),node(1,j)) = 1 - kappa*dt/J(1,j)*(-1/(2*dkesi^2)*(q(2,j)-q(1,j))+1/(2*dkesi*dyita)*(r(2,j)-r(1,j))...
                +1/(2*dkesi*dyita)*(r(1,j+1)-r(1,j-1))-1/(2*dyita^2)*(p(1,j+1)-p(1,j-1)));
    A(node(1,j),node(2,j)) = -2*kappa*dt/J(1,j)*(1/(2*dkesi^2)*(q(2,j)+q(1,j))-1/(2*dkesi*dyita)*(r(1,j+1)+r(1,j)));
    A(node(1,j),node(1,j+1)) = -kappa*dt/J(1,j)*(1/(2*dyita^2)*(p(1,j+1)+p(1,j))-1/(2*dkesi*dyita)*(r(2,j)+r(1,j)));
    A(node(1,j),node(1,j-1)) = -kappa*dt/J(1,j)*(1/(2*dyita^2)*(p(1,j)+p(1,j-1))-1/(2*dkesi*dyita)*(r(1,j)+r(1,j))); 

    A(node(nelx,j),node(nelx,j)) = 1 - kappa*dt/J(nelx,j)*(-1/(2*dkesi^2)*(q(nelx,j)-q(nelx-1,j))+1/(2*dkesi*dyita)*(r(nelx,j)-r(nelx-1,j))...
                +1/(2*dkesi*dyita)*(r(nelx,j+1)-r(nelx,j-1))-1/(2*dyita^2)*(p(nelx,j+1)-p(nelx,j-1)));
    A(node(nelx,j),node(nelx-1,j)) = -2*kappa*dt/J(nelx,j)*(1/(2*dkesi^2)*(q(nelx,j)+q(nelx-1,j))-1/(2*dkesi*dyita)*(r(nelx,j)+r(nelx,j-1)));
    A(node(nelx,j),node(nelx,j+1)) = -kappa*dt/J(nelx,j)*(1/(2*dyita^2)*(p(nelx,j+1)+p(nelx,j))-1/(2*dkesi*dyita)*(r(nelx,j)+r(nelx,j)));
    A(node(nelx,j),node(nelx,j-1)) = -kappa*dt/J(nelx,j)*(1/(2*dyita^2)*(p(nelx,j)+p(nelx,j-1))-1/(2*dkesi*dyita)*(r(nelx,j)+r(nelx-1,j)));
    end
    
    % The following are for corners
    A(node(1,1),node(1,1)) = 1 - kappa*dt/J(1,1)*(-1/(2*dkesi^2)*(q(2,1)-q(1,1))+1/(2*dkesi*dyita)*(r(2,1)-r(1,1))...
                +1/(2*dkesi*dyita)*(r(1,2)-r(1,1))-1/(2*dyita^2)*(p(1,2)-p(1,1)));
    A(node(1,1),node(1,2)) = -2*kappa*dt/J(1,1)*(1/(2*dyita^2)*(p(1,2)+p(1,1))-1/(2*dkesi*dyita)*(r(2,1)+r(1,1)));
    A(node(1,1),node(2,1)) = -2*kappa*dt/J(1,1)*(1/(2*dkesi^2)*(q(2,1)+q(1,1))-1/(2*dkesi*dyita)*(r(1,2)+r(1,1)));

    A(node(nelx,1),node(nelx,1)) = 1 - kappa*dt/J(nelx,1)*(-1/(2*dkesi^2)*(q(nelx,1)-q(nelx-1,1))+1/(2*dkesi*dyita)*(r(nelx,1)-r(nelx-1,1))...
                +1/(2*dkesi*dyita)*(r(nelx,2)-r(nelx,1))-1/(2*dyita^2)*(p(nelx,2)-p(nelx,1)));
    A(node(nelx,1),node(nelx,2)) = -2*kappa*dt/J(nelx,1)*(1/(2*dyita^2)*(p(nelx,2)+p(nelx,1))-1/(2*dkesi*dyita)*(r(nelx,1)+r(nelx,1)));
    A(node(nelx,1),node(nelx-1,1)) = -2*kappa*dt/J(nelx,1)*(1/(2*dkesi^2)*(q(nelx,1)+q(nelx-1,1))-1/(2*dkesi*dyita)*(r(nelx,1)+r(nelx,1)));

    A(node(1,nely),node(1,nely)) = 1 - kappa*dt/J(1,nely)*(-1/(2*dkesi^2)*(q(2,nely)-q(1,nely))+1/(2*dkesi*dyita)*(r(2,nely)-r(1,nely))...
                +1/(2*dkesi*dyita)*(r(1,nely)-r(1,nely-1))-1/(2*dyita^2)*(p(1,nely)-p(1,nely-1)));
    A(node(1,nely),node(1,nely-1)) = -2*kappa*dt/J(1,nely)*(1/(2*dyita^2)*(p(1,nely)+p(1,nely-1))-1/(2*dkesi*dyita)*(r(1,nely)+r(1,nely)));   
    A(node(1,nely),node(2,nely)) = -2*kappa*dt/J(1,nely)*(1/(2*dkesi^2)*(q(2,nely)+q(1,nely))-1/(2*dkesi*dyita)*(r(1,nely)+r(1,nely)));
    
    A(node(nelx,nely),node(nelx,nely)) = 1 - kappa*dt/J(nelx,nely)*(-1/(2*dkesi^2)*(q(nelx,nely)-q(nelx-1,nely))+1/(2*dkesi*dyita)*(r(nelx,nely)-r(nelx-1,nely))...
                +1/(2*dkesi*dyita)*(r(nelx,nely)-r(nelx,nely-1))-1/(2*dyita^2)*(p(nelx,nely)-p(nelx,nely-1)));
    A(node(nelx,nely),node(nelx-1,nely)) = -2*kappa*dt/J(nelx,nely)*(1/(2*dkesi^2)*(q(nelx,nely)+q(nelx-1,nely))-1/(2*dkesi*dyita)*(r(nelx,nely)+r(nelx,nely-1)));
    A(node(nelx,nely),node(nelx,nely-1)) = -2*kappa*dt/J(nelx,nely)*(1/(2*dyita^2)*(p(nelx,nely)+p(nelx,nely-1))-1/(2*dkesi*dyita)*(r(nelx,nely)+r(nelx-1,nely))); 
    
    
    A = sparse(A);
    rhs = phi + dt*mphi;
    phi = A\rhs;
    
    phi = reshape(phi,nelx,nely);
end
    

function s = node(i,j)
    global nelx
    s = (j-1)*nelx+i;
end


